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Abstract 

We study theoretically the complex network of forces that is responsible for 
the static structure and properties of granular materials. We present detailed 
calculations for a model in which the fluctuations in the force distribution 
arise because of variations in the contact angles and the constraints imposed 
by the force balance on each bead of the pile. We compare our results for 
force distribution function for this model, including exact results for certain 
contact angle probability distributions, with numerical simulations of force 
distributions in random sphere packings. This model reproduces many as- 
pects of the force distribution observed both in experiment and in numerical 
simulations of sphere packings. 

Our model is closely related to some that have been studied in the context 
of self-organized criticality. We present evidence that in the force distribution 
context, "critical" power-law force distributions occur only when a parameter 
(hidden in other interpretations) is tuned. Our numerical, mean field, and ex- 
act results all indicate that for almost all contact distributions the distribution 
of forces decays exponentially at large forces. 
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I. INTRODUCTION 



Disordered geometric packings of granular materials have fascinated researchers for 
many years. Such studies, with their applicability to the geometry of glass-forming sys- 
tems, initially were concerned with categorizing the void shapes and densities. More recently, 
partly in recognition of the ubiquity of granular materials and their importance to a wide 
variety of technological processes, interest has focused on how the forces supporting the 
grains are distributed. Visualizations of two-dimensional granular systems demonstrate 
weight concentration into "force chains." It is natural to expect that similar concentrations 
of forces will occur in three dimensions. The distinctive forces in bead packs also give rise 
to distinctive boundary- layer flow and novel sound-propagation properties. |^ 

Ref. 10] presents experiments, simulations, and theory characterizing the inhomogeneous 
forces that occur in stationary three-dimensional bead packs, focussing particularly on the 
relative abundance of forces that are much larger than the average. If the bead pack were 
a perfect lattice, then, at any given depth, no forces would be greater than some definite 
multiple of the average force. At the other extreme, if the network of force-bearing contacts 
were fractal, then fluctuations in the forces (characterized, say, by their variance) would 
become arbitrarily large compared to the average force at a given depth, as the system size 
is increased. Ref. shows that the forces in bead packs are intermediate between these two 
extremes. The forces are unbounded, but the number of large forces falls off exponentially 
with the force. The fluctuations remain roughly the same as the average force, regardless of 
how large the bead pack becomes. A simple model was introduced to understand the results 
of the experiments and simulations. 

This paper presents the detailed analysis of the model introduced in Ref. 0. The model 
yields force distributions which agree quantitatively with those obtained in numerical simu- 
lations of sphere packings. Generic distributions of contacts lead to force distributions which 
decay exponentially at large forces, though a special distribution exists for which the force 
distribution is power law. We discuss the relationship of this model to other related systems 
as well as present the analysis leading to the results that are quoted in Ref. without 
derivation. 

The paper is organized as follows. Section |I| defines the model, discusses several limiting 
cases that have been discussed previously in other contexts, and then presents our analysis of 
the force distribution expected in the context of force chains in bead packs. Special emphasis 
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is placed on one particular contact distribution, the "uniform" distribution, which is the most 
random distribution consistent with the constraint of force balance. We first present a mean 
field solution for this model, and then show that this mean field solution is exact. We also 
obtain exact results for a countable set of non-generic distributions as well as mean-field 
and numerical results for other contact distributions. Evidence is presented that almost all 
contact distributions lead to exponentially decaying force distributions. Section |T| discusses 
numerical simulations of sphere packings, which we analyze to obtain contact probability 
distributions to be used in the g-model. We show that the force distribution predicted by 
the model with this contact distribution agrees quantitatively with the force distribution in 
the simulation. Appendix A presents some mathematical identities concerning the uniform 
g-distribution which are used in the text. 

II. THE q-MODEL 
A. Definition of the model 

Here we introduce the model, which assumes that the dominant physical mechanism 
leading to force chains is the inhomogeneity of the packing causing an unequal distribution 
of the weights on the beads supporting a given grain. Spatial correlations in these fractions 
as well as variations in the coordination numbers of the grains are ignored. We consider 
a regular lattice of sites, each with a particle of mass unity. Each site i in layer D is 
connected to exactly sites j in layer D + 1. Only the vertical components of the forces are 
considered explicitly (it is assumed that the effects of the horizontal forces can be absorbed 
in the random variables qij defined below). A fraction qij of the total weight supported by 
particle i in layer D is transmitted to particle j in layer D + 1. Thus, the weight supported 
by the particle in layer D at the i^^ site, w{D,i), satisfies the stochastic equation: 

w{D + l,j) = l + Y,q,,^{D)wiD,z) . (2.1) 

i 

We take the fractions qij{D) to be random variables, independent except for the constraint 
J^jQij = I5 which enforces the condition of force balance on each particle. We assume that 
the probability of realizing a given assortment of g's at each site i is given by a distribution 
function p^qn, . . . , qny) = {Hj f{<lij)}^(J2j Qij ~ !)• We define the induced distribution 77(g) 
as: 
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= n / dqijpiqu, . . . , Qik = q, ■ ■ ■ , QiN) ■ (2.2) 

Because • • • , Qin) is a probability distribution and J2iLj Qij = 1; the induced distribution 
must satisfy the conditions dqr]{q) = 1, Jq dq q r]{q) = 1/N. 

In this paper we focus on the force distribution Qd{w)i which is the probabihty that a 
site at depth D is subject to vertical force w. We obtain Qoiw) for different distributions 
of g's. We will also consider the force distribution Pd{v) for the normalized weight variable 
V = w/D. For 77(g) = 6{q — 1/N), where each particle distributes the vertical force acting on 
it equally among all its neighbors, the force distribution at a given depth is homogeneous: 
Qd{w) = S{w — D), or Pd{v) = 5{v — 1). At the other extreme, there is a "critical" limit, 
when q can only take on the values 1 or 0, so that weight is transmitted to a single underlying 
particle. For this, as discussed in the next section, the force distribution obeys a scaling 
form and decays as a power law at large forces, Q{w) oc w~'^, where c{N > 3) = 3/2 and 
c{N = 2) = 4/3. We demonstrate that this power law does not occur when q can take on the 
values other than 1 and 0, as is the case for real packings. Generic continuous distributions 
of g's lead to a distribution of weights that, normalized to the mean, is independent of depth 
at large D and which decays exponentially at large weights. We solve the model exactly for 
a countable infinite set of g-distributions, and present mean-field and numerical results for 
other distributions of g's. 



B. The q-model for the "critical" case 

We first consider the case where each particle transmits its weight to exactly one neighbor 
in the layer below, so that the variable q is restricted to taking on only the values and 1. We 
denote this (singular) limiting case of our model by the "go.i limit." Figure shows the paths 
of weight support for a two-dimensional system in this limit. The solid lines correspond to 
bonds for which q = 1. The paths of weight support of particles in the top row are coalescing 
random walks. Since a random walk of length D has typical transverse excursion of D^^'^, for 
the two-dimensional case the maximum weight supported by an individual grain at depth D 
scales as D^^'^. ^ Because D^^'^ ^ D, the mean weight supported at depth it is plausible 
that in the go,i limit the model yields a broad weight distribution. 

The defining equations of the go,i limit of our model are known to be identical to those of 
Scheidegger's model of river networks and a model of aggregation with injection; |PH|JTT[| 



the model is also equivalent to that of the directed Abelian sandpiles. |jT6|,|T5|,|Tj] (The number 
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of neighbors below a particle, A^, corresponds to the dimensionality d in these models.) The 
last equivalence follows [1^,15| if we define Go{Xi] Xq) as the probability that the weight of 
site Xi is supported by site Xq in the same row or below it. The conditional probability 
that Xi is supported by Xq, given that / of the neighboring particles in the row below 
are supported by Xq, is IjN . Thus, 

^ ^ 1 ^ 

Go(Xi; X,) = - Y: G,{X, - el; X^) + 5^^ ^^^ , (2.3) 
i=i 

where {Xi — el} are the neighbors of Xi in the row below it, and the 5-function term follows 
because each particle must support its own weight. Similarly, the probability that two sites 
Xi and X2 in the same row are supported by Xq satisfies: 

X2; Xo) = ^J2J2 GiX, - el, X2 - el; Xo) (2.4) 

i j 

for Xi 7^ X2. These equations are precisely those that describe the behavior of the correla- 



tions of the avalanches in the directed Abelian sandpile. |T6[ |jT9| In this model, an integer 
"height" variable z{X) is assigned each site X on a lattice. The dynamics are defined by the 
rule that if any z{X) exceeds a critical value, Zc, then the variables at m nearest neighbor 
sites along a preferred direction increase by 1, while z{X) decreases by m. In this con- 
text Go{Xi; Xq) is identified with the probability that adding a particle at Xq creates an 
avalanche that topples over the site Xi. Higher order correlations are mapped similarly. 
The distribution of weights in our model is mapped to the distribution of avalanche sizes. 
All these models |pHri|, P^ , P^ have been studied as examples of self-organized criticality. 



12 1 because they lead to power law correlations without an obvious tuning parameter. 



However, in the context of our model, the go,i limit is a singular one, where the probability 
of g 7^ {0) 1} has been tuned to zero. As we shall show in this paper, generic distributions 
77(g), for which the probability that q 7^ {0, 1} is nonzero (no matter how small), yield 
completely different results, with the distribution of weights decaying exponentially at large 
weights. With hindsight, we identify the probability for a river to split in the river network 
model, 1^ and the probability for a colloidal particle to fragment in the aggregation model 
I0| , p!5[] as hidden parameters that were tuned to zero. The corresponding parameter for 



directed Abehan sandpiles is less obvious. 

The equivalence of our model in the go,i limit to the models discussed above can be 
exploited to obtain some results for the distribution of weights. Recalling that the dimen- 
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sionality, d, in these models corresponds to our A^, we know that the weight distribution 
function at a depth D, Qj:){w), has a scaling form for all A^: 

Qniw) = D-'^giw/D") , (2.5) 

where g{x) — > as x — (with a cutoff at w of 0(1)). 

The normalization constraints, dw Qd{w) = 1 and Jl^ dw w Qd{w) = D yield the 
conditions 

a = 6c, 1 + a = 26 , (2.6) 

so that there is only one free exponent. For d = 2, the random walk argument at the 
beginning of this subsection suggests that 6 = 3/2, p which agrees with the exact result. 



11 1 For d > 2, random walks are less likely to coalesce, and this argument breaks down. 



In mean field theory one obtains the analytic result 6 = 2, and exact analytic results 
for directed Abelian sandpiles in all dimensions |]T6| show that mean field theory is valid for 
d > 3, and confirm the result 6 = 3/2 for d = 2. (Our exponent 6 can be identified with 
a + l of Ref. ig.) 

As D ^ (yo, the argument of the scaling function g in Eq. ( |2.5D is small for any finite w. 
Thus, in the go,i limit of our model, the distribution of weights, Q{w), is independent of D 
as — * oo, and is of a power law form, and hence is infinitely broad. 



C. The q-model away from criticality. 

The rest of this paper concerns probability distributions of the g's that do not have the 
property that q takes on only the values 1 and 0. We argue that all such distributions lead to 
force distributions that differ qualitatively from those described in the previous subsection. 
The go,i limit is the only one that yields a power law force distribution; other distributions 
lead to a much faster, typically exponential, decay. In addition, for other g-distributions, 
the distribution for the normalized weight v = w/D, Pd{v) converges to a fixed distribution 
P{y) as D — s> oo. In contrast, in the go,i limit, the quantity Qoiw) converges to a fixed 
function. In this subsection we present evidence for these assertions via both numerical 
simulations and mean field analysis. 
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1. Numerical simulations 



Our numerical investigations all indicate that that for all g-distributions except for the 
go,i limit, the normalized force distribution Pd{v) becomes independent of D as D ^ oo. 
To illustrate typical behavior, we consider the specific g-distribution consisting of — 1 
bonds emanating down from each site with value q = qo < 1/{N — 1) and one bond with 
q = 1 — {N — l)qQ, which has the induced distribution: 

VM = ^Siq - (1 - (iV - l)go)) + ^^(? - Qo) • (2.7) 

Figure displays the normalized force distribution Pd{v) versus v for several different depths 
D in a 3-dimensional fee system (A^ = 3) of dimension 512 x 512 x D, with go = 0.1. Periodic 
boundary conditions are imposed in the transverse directions. As D becomes large, Pd{v) 
converges to a function independent of D which decays faster than a power law. Figure ^ is 
a semilog plot of Pd{v) versus v for several values of showing that the decay of P{v) at 
large D is roughly exponential. To see that this behavior is qualitatively different from that 
of the go,i limit, in Figure |^ we display numerical results for Pd{v) versus v for a system 
which is identical except that go = 0. In contrast to the g = 0.1 case, Pd{v) decays as a 
power law at large v. Also, Pd{v) shows no signs of becoming independent of D as D — oo. 
This is consistent with the result in the previous section that Qd{w) becomes independent 
of D at large D. 



2. Mean field theory 

The technique of the mean field analysis for a general g-distribution is a generalization 
of that used for the go,i case. [|10 



The weight supported by a given site at depth D, Wi{D), depends not only on the weight 
supported by the sites at depth D — 1 but on the values of g for the relevant bonds: 

w,iD + l) = J2q,,w,iD) + l . (2.8) 
j 

In general the values of w at neighboring sites in layer D are not independent; the mean 
field approximation consists of ignoring these correlations. 

As discussed above, when g is allowed to take on values other than and 1, it is useful to 
study the force distribution function as a function of the normalized weight at a given depth, 
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V = w/D. In terms of the normalized weight variable the mean field approximation leads 
to a recursive equation for the weight distribution function Pd{v): 

Pd{v) = Hi / dmi^^) / dv,PD-,{v,)}6{Y^[{D - 1)/D]v,q, - {v - 1/D)) . (2.9) 
i=i ''^ j=i 

The quantity 77(g) is defined in Eq. (|2.2| ). The constraint that the g's emanating downward 
from a site must sum to unity enters only through the definition of 77(g) because there is no 
restriction on the g's for the ancestors of a site. The only approximation here is the neglect 
of possible correlations between the values of v among the ancestors. 

By Laplace transforming, one finds that Pd(s), the Laplace transform of the distribution 
function of the normalized weight Poiv), obeys: 



Pd{s) = e 



-s/D 



/ dqr,{q)Pn-i{sq{D - l)/D) 
Jo 



(2.10) 



Since as -D ^ 00 the distribution Pd{v) becomes independent of D, |21[] Eq. ( |2.10[ ) then 
becomes: 

"1 . -'^ 



P s) 







dqri{q)P{sq) 



(2.11) 



First we show that the weight distribution P{v) decays faster than any power of v for 
all g-distributions except those that only take on the values 0, 1. We expand the Laplace 
transform P{s) in powers of s, P{s) = 1 + Z^^i -Pj-S"'; and plug into Eq. ( p.ll|) , obtaining: 



I + P1S + J2 PjS^ = [1 + ^1^/^" + E PjS^q^)f . (2.12) 

i=2 j=2 

Here, s is the Laplace transform variable, (g-') = Jq dq q^'r](q), and we have used (g) = 
Equating the coefficients of s^ on the left and right hand sides of the equation, we obtain a 
linear equation for Pj: 

P,[N{q^) - 1] = G(P,_i, P,_2, . . . , Pi) , (2.13) 

where G is some complicated polynomial. This can be iterated to obtain Pj for successively 
higher values of j. 

Since Eq. ( |2.13| ) is linear in Pj, Pj can diverge only if its coefficient, [A^(g-') — 1], is zero. 
If g can take on only the values and 1, then (g-') = (g) and [A^(g'') — 1] = for all j > 1- 
However, for any other distribution of g's restricted to the interval [0, 1], the distribution for 
q^ is shifted towards the origin compared to the distribution for g^, whenever j > k. Since 
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{q^) < {q) = 1/N for all j > 2, Eq. ( p.l3| ) has a nonzero coefficient for Pj for all j > 2, which 
means that all moments (v^) of P{v) are finite. (For the special case of j = 1, the equation 
is degenerate; Pi is set by the normalization of v.) If lnP{v) were to behave asymptotically 
as — alnf, then (v^) would diverge for all j > a — 1. Hence P{v) must fall off faster than 
any power of v and {dlnP{v)/d\nv) —oo as f — > oo. 



Now we consider the distribution of weights for non-critical distributions of g's. Moti- 
vated by the geometrical disorder present in granular materials, we focus especially on con- 
tinuous distributions. First we calculate this distribution within a mean field approximation 
for the simplest possible continuous distribution, f{qij) = constant, or pi^qn, . . . ^qi^) = 
{N — l)\S{J2j Qij — 1) (the uniform g-distribution). We show that within mean field theory, 
all "typical" continuous g-distributions lead to a force distribution that decays exponentially 
at large weights. We will show later that the mean field solution is exact for a countable set 
of g-distributions, including the uniform g-distribution. 



One example of a g-distribution that can lead to an exponentially decaying distribution 
of weights is the "uniform" distribution of g's, for which the probability of obtaining the 
values qn, . . . , qiN is p(gii, • • • , Qin) = {N — lij — !)• We show in Appendix A that 

this distribution induces riu{q) = {N — 1)(1 — g)^^^. Thus, for this g-distribution in the limit 
D ^ oo the mean field force distribution is the solution to the self-consistent equation: 



D. Weight distributions away from criticality: Mean field results 



1. Mean field theory for the uniform distribution 




(2.14) 



First consider N = 2. For this case ^^(g) = 1, so Eq. ( ^.141 ) becomes: 




(2.15) 



Letting V{s) = (P(s)) 



and u = qs, one obtains: 




(2.16) 



Differentiating with respect to s yields: 
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as 



(2.17) 



which can be integrated to yield 



Vis) 



l~Cs 



(2.18) 



The constant of integration C is determined by the definition of the mean, /q°° dv v P{v) = 
= 1. Thus, C = %\s=o = -1/2. Hence one finds P{s) = 4/(s + 2f and P{v) = 

This method can be generahzed for all A^. Defining Vn{s) = {Pn{s)Y^^ , inserting in 
Eq. ( p.l4| ), and differentiating — 1 times, one finds that Vn{s) obeys the differential 
equation: 



ds^- 



-(.--V^(.)) = (iV-l)!V^-(.) 



(2.19) 



A solution to this equation is Vn^s) = where C is any constant. This can be shown by 
induction: Assume that 

c ^^-1 



jN-2 / 



s + C 



{N-2)\ 



s + C 



(2.20) 



Then 



d 



N-l 



ds^- 



1 c 



s + C J ds 



jN-1 , fi 

^-^A{s + C-C)s^--^ 



-C 



ds 



{N -2)\ 



C 



s + C 

c 



S + C 



(2.21) 
(2.22) 

(2.23) 



.s + C, 

Since direct substitution can be used to show that the identity holds for A^ = 2, it holds for 
all N. 

The condition ^ = is satisfied when C = N. Hence one finds the weight distribu- 
tion: 



Pn{v) 



(N-iy. 

The question of uniqueness of this solution is discussed below. 



(2.24) 
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2. Mean field asymptotic force distribution for generic continuous q- distributions 



We now show that, within mean field theory, generic continuous g-distributions lead to 
weight distribution functions P{v) for the normalized weight v which have the asymptotic 
forms P{v) oc v^^^e''^^ as v ^ oo and P{v) oc v^~^ as v 0. 

We consider g-distributions of the form p^qn, . . . ,qiM) = {Hj f{(lij)}^(J2j Qij — 1) (the 
uniform distribution is f{qij) = constant). If f{qij) has a nonzero limit as qij 0, and 
does not have a 5-function contribution at q^j = 0, then phase space restrictions imply that 
the induced distribution 77(g) ~ (1 — g)^~^ for g — 1. This is because if a site receives a 
fraction q of the weight from one of its predecessors, then the fractions received by all the 
other successors of that predecessor, {g2 . . . gAr} must add up to 1 — g. For g close to 1, this 
gives a phase-space volume of the order of (1 — g)^~^. 

To determine the large v asymptotics of P{v), we use the result of Sec. ( [11 C 2| ) that P{v) 
must fall off faster than any power of v. We write the D ^ 00 limit of Eq. (|2.9|) as: 



P(v) 



N 

n 



[ / dv,F{v,) 5{v - E 

1^0 J j 



^0) 



(2.25) 



dq,jP{vj/qj)ri{qj)/qj . 



(2.26) 



Since P{v) decays quickly (in particular, faster than 1/f), the apparent singularity near 
g = in Eq. (p.26|) is not really there. The integral is dominated by g ~ 1. This follows 
because 

(9 In P(t>/g) 



P{v/q) = P{v) exp 

= P{v) exp 
= P{v) exp 



dq 

V din P{v/q) 



q=l 



q^ d{v/q) 
dliiPiu] 



d\n.u 



1-g) 



(2.27) 



Since dhi P{u) /dlnu —00 asu 00, this expression becomes very small as 1— g increases. 
Thus, for large v, since 77(g) ~ (1 — q)^~^ for g ~ 1, 

N-l 



F{v) ~ P{v)/ 



dlnPiv] 



dlnv 



(2.28) 



Already it is clear that P{v) for any generic g-distribution has the same large-v asymp- 
totics as the uniform distribution, since the asymptotics are determined entirely by the phase 



11 



space restrictions on 77(g) for g ^ 1. This decay also can be demonstrated explicitly by as- 
suming faster and slower decays and showing inconsistency with Eq. ( p.25| ). If P{v) were 
to decay faster than exponentially, then the convolution in Eq. ( |2.25| ) would be dominated 
by the region where all the v'jS are roughly equal. But since [P{v/N)]^ ^ Eq. ( p.25| ) 

cannot be satisfied. On the other hand, if P{v) were to decay slower than exponentially, 
then the convolution would be dominated by the region where one of the fj's is ~ f and the 
others are 0(1). Eq. (|2.25| ) would then imply: 

N-l 



P{v) ~ P{v)/ 



d\nP(v) 



d\nv 



(2.29) 



Since the expression in square brackets diverges with v, this is not possible either. Thus 
one must have P{v) = h{v) exp[—av], where h{v) varies more slowly than an exponential. 
Eq. ( |2.25| ) then implies: 



hiv) 

This is satisfied by h{v) 



^ 1 I 



j'^'^, so that 



(2.30) 



P(v) 



V 



N-l 



exp[— Of] 



(2.31) 



for V 00. 

Hence we have shown that for generic continuous g-distributions, within mean field theory 



P{v) 



,N-l 



exp(— at>) Bs V 00. 



3. Mean field theory for singular q- distributions. 

We have shown that all g-distributions which satisfy the condition f^dq 77(g) ~ (1 — 
g)^~^ as g 1 have a weight distribution within mean field theory that is of the form 
P{v) ~ f^~^exp[— Of] for large v. This condition on 77(g) is satisfied under fairly general 
assumptions: one requires (1) that the probability density for any qij in Eq. ( |2.8| ) have 
a nonzero g^j limit and (2) that it not have a ^-function contribution at g^j = 0. 
However, as we shall see below, to compare the results of the g-model to molecular dynamics 
simulations and to experiments on real bead packs, it is useful to consider the case where 
there is a finite probability for some of the g^'s to be zero, which implies that the induced 
distribution 7/(g) has a 5-function at g = (and in some cases at g = 1). |^ Such a choice 
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for ri{q) is also useful in examining the crossover from the critical go,i limit to the smooth 
g-distributions considered in the previous subsection. We will see that g-distributions of this 
type lead to force distributions P{v) that decay exponentially, though with different power 
laws multiplying the exponential than for continuous g-distributions. 

We first note that, when ri{q) has a finite weight at g = 1, it is impossible for P{s) to 
diverge at any s. The solutions of the form P{s) ~ l/(g + s/ sq)^ obtained in section ( |11 D 1| ) 
were possible because, in Eq. ( p.llD , the integral over g reduces the singularity, which is 
compensated by the exponentiation. With a finite weight at g = 1, close to a divergence at 
So one would have P{s) oc [P{s)]^ , which would be impossible as s ^ sq. 

It is instructive to consider first a simplified version of such singular g- distributions. Let 
us consider the case of = 2, and assume that r]{q) has the form: 

r^(g) = 1(1 - 0){6{q) + 6(1 - g)} + e6{q - 1) , (2.32) 

with < 6* < 1. This ri{q) satisfies the conditions j dq rj{q) = 1 and J dq q ri{q) = 1 /2 for all 
6. Eq. (|2.11|) then simplifies to: 



Pis) 



l{i-e) + ^{i-0)p{s) + ep{s/2)' 



(2.33) 



where we have used the fact that -P(O) = 1. Eq. (|2.33|) can be solved as follows: for small 
s, we know that P{s) = 1 — s + O(s^) (the coefficient of the linear term being fixed by the 
normalization condition J dv v P{v) = 1). Starting with a small negative value of s, where 
P{s) is approximated as 1 — s, Eq. (|2.33| ) can be iterated to find P(2"s) for n = 1, 2, . . . 
(the correct root of the quadratic equation is chosen by requiring P{s) = 1 for P{s/2) = 1). 
Eventually the result of this iteration scheme is complex rather than real, signifying that 
s is in a region where P{s) has a branch cut. It is easiest to find the origin sq of this 
branch cut by adjusting P{sq/2) so that Eq. ( |2.33D has a double root for -P(so), and then 
iterating backwards to obtain P(so/2"). As n — > oo, by matching on to the requirement that 
P{s) = 1 — s for s ^ 0, one can obtain sq. It is clear from Eq. ( ^.33] ) that in the vicinity of 



So P{s) is of the form -P(so) + ay/s — sq. This yields 

P{v) ~ v'''^^^ exp[-sov] for f ^ oo . (2.34) 

Although the power law prefactor is different from that in Eq. ( p.31| ), there is still an 
exponential decay. 

We now consider possible changes to Eq. ( 2.34 ) from choosing //(g) of a more complicated 
form than Eq. ( p. 321) . For any T]{q) of the form 
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riiq) = E - + (1 - E c,)6{q)6{q) (2.35) 



i=0 



with < A < 1, one can use the method outhned above to find that P{s) has a square-root 
branch cut at some Sq. This answer is not affected by making n large, so long as Cq remains 
nonzero. As n — oo, with all the q's for i > tending to zero, we can approach arbitrary 
continuous distributions for ri{q) with 6 functions at g = and q = 1- 

For N > 2, Eq. ( p.33| ) is changed to a higher order equation. This, however, does not 
generically change the results above. Even for higher order equations, the degeneracy of the 
roots generally occurs only pairwise, so that close to the point of degeneracy the singularly 
ranging roots still have a square-root singularity. It will, however, be possible to find non- 
generic choices for r](q) that could result in an asymptotic form P{v) ~ v^^^^^^ exp[—av] 
with N > m> 2. 



E. Beyond mean field theory 

1. Proof that mean field theory is exact 

In this section we prove that the mean field solution presented in the previous subsection 
is an exact solution of the model with the uniform g-distribution for any A^. 

In general, the mean field theory presented above is not exact because it does not account 
for the fact that two neighboring sites in row D + 1 both derive a fraction of their weight from 
the same site in row D. Suppose a site j in row D + 1 has w{j) much larger than the average 
value. Then it is likely that the weight supported by an ancestor w{i) in row D is larger than 
average also. Because this ancestor transmits its weight to a neighboring site in row D + 1 as 
well, there is a "correlation" effect that creates a greater likelihood that in a given layer sites 
supporting large weight are close together. On the other hand, there is a "anticorrelation" 
effect arising because J2j Qij = 1; if a large fraction of the weight from site i is transmitted 
to site j, then small fractions are transmitted to the other "offspring" sites. When the g's 
are chosen from the uniform distribution, these "correlation" and "anticorrelation" effects 
cancel exactly. 

The result that the mean field correlation functions are exact for the uniform distribution 
of g's can be understood by considering the system in terms of weights on bonds. Each 
bond {ij} corresponds to a particle with "energy" Eij = Viqij. Moving down by one layer 
corresponds to having groups of particles colliding at each site and emerging with different 
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energies, subject to the constraint that the total energy of all particles colliding at each 
site is unchanged by the collision. 

For the "uniform" g-distribution, each collision takes particles of energies Ca^, . . . , Cq,^ 
and changes their energies to Ea^, . . . , -Eq,^, subject only to the constraint that J2^a = J2 Ea- 
If we start with a "microcanonical" ansatz for the phase-space density, i.e. that it is uniform 
over the space J2 Ea = E, then it is preserved by the collisions. Hence, the microcanonical 
density is the correct one for this system. 

With a microcanonical density for a large collection of particles, the density for any finite 
subgroup is canonical (in the thermodynamic limit). p6| Thus, we have shown for this case 
that the distribution of "bond forces" is exponential, which is the most random distribution 
consistent with the constraint that the sum of the forces is fixed. |2^j27[ 



Note that this argument does not hold for g-distributions other than the uniform one. 
For instance, in the go,i limit, each collision takes all the energy of the group and gives it to 
one of the colliding particles. Thus, even if we start with the microcanonical distribution, 
it breaks down at the very first step. For general g-distributions, the phase space density is 
not separable, i.e., mean field theory is not exact and there are spatial correlations within 
each layer. 

The explicit algebraic proof proceeds by constructing exact recursion relations for the 
correlation functions describing the weight distribution in the model in row D + 1 in terms 
those for row D, and showing that the mean field correlation functions are invariant under 
this recursion. We ignore the weight added in each row because we are looking for the fixed 
distribution very far down the pile. 

Let Poiui) be the probability that site i in row D supports weight Wj, Pniui^^Ui^) 
be the probability that sites ii and 12 support weight and Ui^, respectively, and 
Pj:,{uij^,Ui2, . . . , Ui^) be the normalized joint distribution describing the probability that sites 
ii, 12, . . . ,in support weights Ui^^, . . . ,Ui^, respectively. The mean field joint probability dis- 
tributions are given by the mean field P{u) and 

P(mi, ^2, ■ ■ ■ , Un) = P{u^)P{u2) ■ ■ ■ P{Un) . (2.36) 

Consider the M-point correlation function in row D + 1 that is obtained when all the 
correlation functions in row D are the the mean field ones. Let {ui} be the weights in row 
D and {vi} be the weights in row D + 1. Consider a cluster of sites j = 1, . . . , M in row 
D + 1, with ancestors in row D at sites i = 1, . . . ,p. (The labels do not imply any particular 
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spatial relation of the sites.) The g's describing the bonds emanating from ancestor i are 
qu, where / = 1, . . . , A^. We define rju^j) to be 1 if sites i and j are connected by bond il 
and zero otherwise. The M-point correlation function in row -D + 1, Pd+i{vi, . . . , vm), must 
obey: 

N 



Pd+i{vi,. . . ,vm) = T\\ dqa... dqiN {N - 1)\ 6{1 - J^Qik) 
i=l ^ ^0 ^0 

roo M / V N \ 

/ duiPoiui)^ \{5{vj ■ (2.37) 

■'^ j=i \ 1=1 1=1 J 

We define the general Laplace transform 

/■OO nOO 

P{SI, ...,Sn)= dVi--- dVnP{Vu • • • , Vn) 6- . (2.38) 

Jo Jo 
Laplace transforming Eq. (|2.37|) , one obtains: 



Pd+i(si, . . . , Sm] 



p n rl ^ 

n / dqa ... dq,M {N - 1)\ 5(1 - ^ q,k) 

^=l ^0 ^0 fc^i 



M N 



PD\^Y>a{j)quSjj . (2.39) 

For Pd{x) = (1 + x/N)~^ , one can use the condition Y^fLi Qu = 1 to write: 

Pd+i(si, . . . ,sm) = n / dqn... dq^ (N - 1)\ 6{1 - J2 lik) 
1=1-'^ "'0 k=i 

N M \ 

5:g,Kl + E^'Kj>,/iV)l • (2.40) 

^1=1 j=i 

Using the identity 

n (a„)-^ = {N-1)\ f'dx,... [' rfx^iil^^l^-ll^M^ (2.41) 
Jo Jo [aiXi + . . . + aNXN) 

with a„ = 1 + Y.jLiVin{j)sj/N, one finds: 

p ^ I 

If a given bond {in} connects to no sites in the descendant cluster, then every term in the 
sum in the denominator of Eq. ( p.42| ) is zero, and the {inY^ term in the product is unity. If 



the bond connects to a site in the descendant cluster, then riin{j) is unity for exactly one j. 
Each site j in the descendant cluster is connected to exactly antecedents in row D, so: 

16 



M 



Po»(s^,...,SM) = nj^^^. (2.43) 

Thus, the mean field correlation functions are preserved from row to row for this q- 
distribution. 

2. Other q distributions 

We have identified a countable set of g-distributions for which mean field theory is exact, 
those of the form /(%) = q^, for all integer r (the uniform distribution is r = 0). The result- 
ing force distribution Pr{v) oc has Laplace transform Pr{s) = (1 + s/{Nr))~'^^. 
The demonstration that this solution is exact follows precisely the same line of reasoning as 
for the r = case presented in the previous subsection, utilizing the identity: p4 



11 K = rp. xijv / ^^1^1 •••/ dXNX^ . 2.44 

n=i W\T)\ -^0 -^0 (ai^i + • • • + aNXNy 

In terms of the particle collision picture discussed at the beginning of Sec. ( [II E 1|) , a 



general value of r corresponds to the particles having an energy which is the sum of r + 1 
components (which may be viewed as as spatial coordinates in some underlying space), 
each one of which is conserved individually in a collision. The microcanonical phase-space 
density, uniform in the N{r + l)-dimensional space, is preserved by the collisions, and yields 
the Pr{v) state here. 

The result that mean field theory yields an exact solution of the model holds only for a 
very limited class of g-distributions. For general g-distributions, the phase space density is 
not separable, i.e., mean field theory is not exact and there are spatial correlations within 
each layer. For example. Figure § shows P{v) for a two-dimensional system with N = 2 
and the g-distribution where the two bonds emanating from each take on the values go 
and 1 — go, with go = 0.1. In the model, a site {i,D) is connected to sites {i,D + 1) 
and {i + 1 mod L,D + 1); in the mean field calculation, site {i,d) is connected to sites 
{pi{i),D + 1) and {p2{i),D + 1), where pi and p2 are permutations of (1,...,L). This 
method of simulating mean field theory destroys the spatial correlations between ancestor 
sites, while ensuring that every site has exactly two ancestors and two descendants. The 
numerical data were obtained by averaging P{v) for rows 10, 001 — 20, 000 in a system of 
transverse extent L = 20, 000. This figure demonstrates explicitly that the mean field force 
distribution P(v) is not exact for this distribution. However, the deviations of the mean 
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field theory from tlie direct simulation are extremely small for v ^0.1, so mean field theory 
provides an accurate quantitative estimate for P{v) over a large range of v. 



3. Uniqueness of the steady state distribution 

In this subsection we show that our results (numerical and analytical) for the force 
distribution do not depend on either the boundary conditions imposed at the top of the 
system or on the specific realization of randomness a particular system might have. 

Consider a system of finite transverse extent L, in which weights {w{i)} are put on the 
particles in the top row. The weight then propagates downwards according to Eq. ( |2.8| ). If 
we now consider the same system with a different loading on the top row, {w{i) + 5w{i)}, 
then since Eq. ( |2.8| ) is linear in w, the difference between the two solutions satisfies the 



homogeneous equation 

6w{D + = Y.q,^,{D)6w{D,i) . (2.45) 

i 

Summing up both sides, we see that J2j Sw{D + 1, j) = J2i Sw{D, i), which means that the 
total excess weight placed on the top of the system progagates downwards unaltered. Such 
a change only affects the normalization of our distributions. Thus, if we are interested in 
the normalized distribution Pci(f), we can without loss of generality consider perturbations 
{6w{D,i)} satisfying the constraint J2i^'w{D,i) = 0. 

Equation ( |2.45|) can be viewed as a two stage process: (1) each Sw{D,i) splits into 
parts, qij{D)6w{D,i), and (2) the fragments qij{D)6w{D,i), with i running over the 
neighbors of j in the row above it, combine to give 6w{D + 1, j) . The important thing is that 
all the qi,/s are positive. Thus if we define the total difference between the two configurations 
as A{D) = J2i \Sw{D, then because all the g's are positive and J2j lij = 1, A is unchanged 
in the first step, while in the second step it can either stay constant or decrease (depending 
on whether the signs of the fragments are the same or different). Further, while for any 
particular value of D it is possible for A(D) to be equal to A{D + 1), the only way in which 
A(D) can remain unchanged as D increases is if all the positive Sw^s are segregated from all 
the negative ones. Even if this is the case in the top row, this becomes increasingly unlikely 
as D is increased. In fact, if the minimum distance between positive and negative 6w is / in 
the top row, and if there are no 5-functions in rjlq), then A{D + 1) must be less than A(D) 
for D > I. Thus for a system of finite transverse extent, the distribution of weights at the 
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bottom of the system is independent of the loading on the top row in the hmit that the height 
of the system is infinite. For the case when all dimensions of the system are made infinite 
the situation is trickier; due to the conservation of J2 under the evolution of Eq. ( p.45| ) 
discussed above, if one were to make 5w positive on one side of the top row and negative on 
the other half, for a system of transverse extent L it would require a height 0{L'^) for the 
effects of this perturbation to "diffuse" away. For generic loading at the top, however, we do 
not expect such an anomalous concentration of fluctuations into only the longest wavelength 
modes of the system, and A(D) should decay with D even if all dimensions of the system 
are enlarged. 

We have seen that the distribution of weight at the bottom of any infinite system is 
independent of the details of how forces are distributed at the top, at least in the limit 
when the height of the system is taken to infinity before its transverse dimensions. This is 
true for each system individually, and is therefore also true for the full ensemble of systems 
with different realizations of randomness (the choice of g^j's), so that the solutions we have 
obtained so far for quantitites like P{v) are unique. For any particular system, however, 
the weights on the different sites at the bottom do depend on the gjj's; in fact, with all 
the gjj's specified, the weights on the different sites are completely determined. Even for 
a single system, however, statistics can be obtained by measuring quantities across all the 
sites in the bottom row; for a system of infinite transverse size the measurements then lead 
to distributions. At least for the "uniform" g-distribution, any quantity like, say, P{v), is 
the same, whether obtained by averaging over sites in a single system or for a single site over 
the entire ensemble. This is because, as we have seen, the ensemble averaged distribution of 
{f 1, . . . , vl} across the bottom row is of the form P(f i, . . . , vl) = P{yi)P{y2) ■ ■ ■ P{vl). For 
any single system chosen randomly from the ensemble, this is the probability density that 
the normalized weights in the bottom row take on the specific values {fi,f2, . . . The 
probability that / of these L sites will have Vi greater than some Vq is then 



as L oo, Z/L is sharply peaked around P{v)dv, so that the site averaged result is 
the same as the ensemble average. We expect this to be the case even for more general 
g-distributions, for which the ensemble averaged P{vi,V2, ■ ■ ■ ,vl) does not have a product 
form, so long as the transverse correlation lengths are finite. 




(2.46) 
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III. NUMERICAL SIMULATIONS OF SPHERE PACKINGS 



We now discuss the relevance of the g-model to granular materials. Although we have 
shown that the g-model yields an exponentially-decaying force distribution independent of 
the details of the g-distribution, to make quantitative comparison of this model to granular 
systems, we must know the g-distribution for a granular material. To make this comparison, 
we have performed molecular dynamics simulations of three-dimensional sphere packings, 
analyzed the contact distributions to estimate the distribution of g's, and then calculated 
the force distribution in the sphere packing and compared it to that predicted by the g- 
model. Our simulations yield results for the contact force distributions that are consistent 
with previous work; p5|-pT|] the new ingredient here is that the geometry of the packing is 
characterized simultaneously, allowing testing of the statistical assumptions underlying the 
g-model. 

Our simulation consists of 500 spherical beads of weight and diameter unity in a uniform 
gravitational field with gravitational constant g = 1, interacting via a central force F of the 
Hertzian form, F = Fo(5r)^/^. Here, Fg is the force constant, chosen so that a sphere has 
a deformation of 6r = 0.001 when subjected to its own weight, and Sr is the deformation 
of each bead at the contact. The box containing the beads had a fixed bottom, and lateral 
dimensions of 5.5 x 5.5. In each simulation, the spheres are initially placed in a loose 
rectangular lattice with lattice constants of 1 x 1 x 1.5 and have random initial velocities 
uniformly distributed in the range —Vmax < Vx,y,z < Vmax, where Vmax = 50 is large enough 
to yield significantly different packings from run to run. By freezing the motion of the 
beads whenever the total kinetic energy of the system reaches a maximum, the kinetic 
energy of the system is reduced and eventually the spheres all settle to the bottom of the 
box. Starting with a flat bottom, the regularity of layer-like packing reduces as the height 
increases. A rough bottom was obtained by selecting the beads with height between H and 
H + 1 (typically, H ~ 10) and this rough bottom was used for the next simulation. Within 
a few iterations, the statistical properties of the rough bottom becomes independent of its 
initial configuration; this configuration of spheres at the bottom of the box is then fixed and 
used as a boundary condition for subsequent packing simulations. 

In our packings, a sphere can have up to 6 contacts on its bottom half. However, on 
average, the three strongest vertical forces at these contacts sustain over 98% of the load; 
three or fewer particles supported at least 90% of the weight for over 92% of the particles. 
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Therefore, comparison with the g-model with = 3 is reasonable. 

We estimate the g-distribution for the sphere simulation by calculating the fractions of 
the total vertical force supported by each of the three strongest contacts. To display 
our results for the g-distribution for the simulation of hard spheres, we define the variables 
c^i = {Qs ~ Q2)/V^, 0L2 = qi- [^] Because Y^^^iQj = 1; the possible values of the a's can 
all be represented as points in the interior of an equilateral triangle, where the values of 
q are the perpendicular distances to each side of the triangle. Moreover, for the uniform 
g-distribution, the density of points in the triangle is constant. If one orders the g's so that 

> (I2 > Qi, then in terms of the a variables, all the points lie in the triangle shown in 
Figure |^, which is bounded by the lines ai > 0, 02 > 0, and ^/Sai + 0:2 < 1. As Figure ^ 
demonstrates, there is some deviation from the uniform g-distribution because a nonzero 
fraction of the particles have gi = 0:2 = 0. A reasonable description of the numerically 
observed particle contact distribution is obtained by taking each particle and assigning with 
probability p, I, and u into "point" , "line" , and "uniform" pieces. In the "point" piece one 
of the g's has value unity, and the other two are zero. In the "line" piece one of the g's is 
set to zero, and the other two are determined as in the N = 2 uniform distribution. Finally, 
the particles in the "uniform" piece have their g's determined exactly as in the A^ = 3 
uniform distribution. Our numerical data for the spheres are consistent with the values 
p = 0.017 ± 0.0023, / = 0.1635 ± 0.007, u=l-l~p = 0.8195 ± 0.007. 

We now discuss our results for the distribution of vertical forces. First, we calculated 
the force distribution at several different depths D. Our numerical data indicate that if 
one considers the normalized force v = w/D, the force distribution P{v) indeed becomes 
independent of depth for D ^ 5, and it decays exponentially at large v. The data were 
obtained by making a histogram of the vertical force exerted by spheres in horizontal slices 
of width AD = 1. The scales are set by the normalization requirements /q°° dv P{v) = 1, 
J,^dvvP{v) = l. 

We now compare the results of these molecular dynamics simulations to results from the 
g-model. Figure |^ shows P{v) calculated via numerical simulation of the g-model Eq. ( p.l| ) 
with /(g) = 1 at depth D = 1024 on a periodically continued fee latice of side 1024, with 
A^ = 3. Within our approximation of placing the grains on a uniform lattice, the reasonable 
choice for A^ is the dimensionality d of the system: For d = 3 the grains are approximated 
as being in triangular lattice layers, with each layer staggered relative to the next, so that 
each grain has three neighbors. As expected, since the mean-field distribution is exact for 



21 



the uniform case, there is excellent agreement with Eq. (|2.24| ). On the same graph we show 
P{v) obtained in the sphere simulation described above. Both the sphere simulation and the 
g-model exhibit a P{v) that decays exponentially at large v. The quantitative agreement 
between the two is surprisingly good considering the "arching" |^ in the sphere simulation, 
as reflected in the "line" and "point" pieces of the g-distribution for the spheres. To examine 
the effects of arching on the results, we examined the force distribution resulting from the 
"g-model" with the three-piece g-distribution, which more closely approximates that of the 
sphere simulation. Figure ^ shows the numerically calculated P{v) for the g-model with 
the three-piece distribution with p = 0.017 and / = 0.1635, together with the solution for 
the uniform distribution and the numerical data from the sphere simulation. Changing the 
g-distribution has little effect on P{v); to the extent that there is a change, it appears to 
improve the already good agreement between the g-model and the sphere simulation. 

Thus, our simulations indicate that our sphere packings are reasonably well-described 
(at the ~ 15% level) by the uniform g-distribution. Deviations from this g-distribution are 
observed; accounting for them improves the already good agreement between the g-model 
and the simulations. 



IV. DISCUSSION 



This paper presents a statistical model for the force inhomogeneities in static bead packs 
and compares the results to numerical simulations of disordered sphere packings. The ir- 
regularities of the packing are described probabilistically, in terms of spatially uncorrelated 
random variables. Although there is a special g-distribution for the g-model that leads to 
a force distribution that decays as a power law at large forces, we have presented evidence 
that the force distribution decays exponentially at large forces for almost all g-distributions. 
We obtain exact results for all the multipoint force correlation functions at a given depth 
for a countable set of g-distributions, including one that is "generic" (the "uniform" dis- 
tribution). The force distribution function for the uniform case agrees quantitatively with 
that obtained for the sphere simulation. Our numerical calculations demonstrate that a 
modifled distribution of g's which more closely approximates that observed for the sphere 
simulation improves the already good agreement between the force distribution predicted by 
the g-model using the uniform g-distribution and simulations of spheres. Thus, this model 
appears to contain some essential features of the force inhomogeneities in granular solids. 
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Neither our simulations nor the q-model of Eq. ( |2.1|) captures all features of real bead 
packs. In our simulations, we have included only central forces and have ignored friction; 
the q-model ignores the vector nature of the forces, assuming that only the component along 
the direction of gravity plays a vital role. The qualitative consistency between the results 
obtained using the different methods as well as with experiment [0 provides some indication 
that the effects that we have neglected do not determine the main qualitative features of 
the force distribution at large v. 

Several avenues for future investigations are evident. It should be straightforward to 
extend the analysis of the model to calculate longitudinal (along the direction of gravity) 
correlations of the forces. It is not obvious how to measure these correlations experimentally, 
but comparison to sphere simulations is clearly possible and would provide further tests of 
the statistical model. Similarly, the theory makes clearcut predictions for the multipoint 
correlation functions, which can be tested both by experiment and by simulations. The 
model can be generalized to apply to a broader variety of situations by including vector 
forces as well as incorporating boundary effects. Most interestingly, we plan to investigate 
whether the statistical theory developed here can be extended to provide new insight into 
the complex dynamical effects exhibited by granular materials. [|I| 

In summary, we have presented and analyzed a statistical model for force inhomogeneities 
in stationary bead packs. The model, which predicts that force inhomogeneities decay 
exponentially at large forces for almost all contact distributions, agrees well with numerical 
simulations of sphere packings as well as experiment. 
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APPENDIX A: THE UNIFORM q-DISTRIBUTION 



Here we consider the "uniform" g-distribution, which is the simplest g-distribution consis- 
tent with the restriction that J2i=i Qi — 1- It is obtained by choosing each of gi, 52, ■ ■ ■ , Qn-i 
independently from a uniform distribution between and 1, setting — 1 — J2f=i Qj, and 
then keeping only those sets where qn > 0. Here we show that rju{q) — j^i^ — q)^"'^ for 
this distribution. 

For N = 2, a one chooses qi between and 1, then q2 = I — qi must also be between 
and 1, so that i]u{q) = 1. When = 3, configuration will be retained only if q'l + 52 + ?3 ^ 1- 
Therefore, the probability of obtaining a value of q is given by: 

rjuiq) =M [' dqi /' dq2S{l - q^ - q^ - q) = M " dqi = M(l - q) , (Al) 
Jo Jo Jo 

where M is a normalization constant. Since dqr]{q) — 1, one immediately finds M — 2 — 
(N-l). 

For general N, riu{q) can be written: 

Vu{q)^V{q)/ rV{q)dq, (A2) 
Jo 



where 

dqnS{l -q- 

1=2 

Using the identity 

/ . Qm) dQn-k+l — T 
m=l ^ m=l 

one can show that 

\N-2 



V{q) = dq2 dqs... dqj{l - g - E ^0 ■ (A3) 

Jo Jo Jo 



~l_y^" n—k+l 1 n—k 



m(q)^(N-l){l-qr-\ (A5) 
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FIGURES 



FIG. 1. Schematic diagram showing the paths of weight support for a two-dimensional system 
in the go,i limit where each site transmits its weight to exactly one neighbor below. The numbers 
at each site are the values of w{i, D). 



FIG. 2. Linear- linear and log-log plots of the normalized weight distribution function Pd{v) 
versus v for a three-dimensional system on an fee lattice {N = 3), for the g-distribution defined 
in Eq. ( p.7[ ) with go = 0.1. The distribution Pd{v) appears to become independent oi D as D 
becomes large, and decays faster than a power law at large v. 



FIG. 3. Semilog plot of the normalized weight distribution function Pd{v) versus v for a 
three-dimensional system on an fee lattice (A^ = 3), for the g-distribution defined in Eq. ( |2.7D 
with go = 0.1. The behavior of Pd{v) at large v is consistent with exponential decay. 



FIG. 4. Linear- linear and log-log plots of the normalized weight distribution function Pd{v) 
versus v for a three-dimensional system on an fee lattice (A^ = 3), for the g-distribution defined in 
Eq. ( |2.7| ) with go = 0. For this special case, the distribution Pd{v) does not become independent 
oi D as D becomes large. The asymptotic decay of Pd{v) at large v is a. power law. 



FIG. 5. Comparison of force distribution P{v) versus v for simulation of mean field theory and 
of the original model equations ( p.l[) for a system with the g-distribution Eq. ( ^.T]) with N = 2 and 
go = 0.1. Both data sets were obtained by averaging the bottom 10000 rows of a 20000 x 20000 
system. Mean field theory does not yield the exact P{v) for this g-distribution. Nonetheless, it 
provides an accurate quantitative estimate for P{v) over a broad range of v. 
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FIG. 6. Scatter plot of contact variables ai, 02 (defined in the text) obtained from the 
sphere simulation described in the text. The graph has 3229 points. On this plot, the uniform 
(/-distribution would have a uniform density of points. The "arching" in the simulation is reflected 
in the fact that a nonzero fraction of points have a2 = 0. 



FIG. 7. The distribution of forces P{v) as a function of normalized weight v = w/D at a given 
depth D. Dashed line: P{v) at D = 1024 obtained via numerical simulation of the model Eq. (2.1) 
with f(q) = 1 on a periodically continued fee lattice of transverse extent 1024. Solid line: Puiv) 
obtained from the analytic mean field solution, Eq. ( |2.24| ). The points are P{v) obtained in the 
sphere simulation described in the text at depth D = 10 (triangles) and averaged over depths 
D = 6 through D = 13 (diamonds). There are no adjustable parameters; the scales are set by the 
normalization requirements dv P{v) = 1, /q°° dv v P{v) = 1 



FIG. 8. The distribution of forces P{v) as a function of normalized weight v = w/D at a given 
depth D. Dashed line: Puiv) obtained from the analytic mean field solution for q-model with 

= 3, Eq. ( 2.24 ). Solid circles: P{v) obtained in the sphere simulation averaged over depths 
D = 6 through D = 13. Open squares: P{v) at -D = 16 obtained via numerical simulation of the 
model Eq. (2T) with the three-part g-distribution described in the text, with parameter values 
given on the graph, on a periodically continued fee lattice of transverse extent 256. This figure 
demonstrates that using the measured g-distribution instead of the uniform g-distribution improves 
the already good agreement between the g-model and the sphere simulation. 
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